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Abstract 

Discretization methods for ordinary differential equations based on the use of matrix expo- 
nentials have been known for decades. This set of ideas has come off age and acquired greater 
urgency recently, within the context of geometric integration and discretization methods on 
manifolds based on the use of Lie-group actions. 

In the present paper we study the approximation of the matrix exponential in a particular 
context: given a Lie group G and its Lie algebra g, we seek approximants F{tB) oi eicpitB) such 
that F{tB) £ G if -B G g. Having fixed a basis Vi, . . . , Vd of g, we write F{tB) as a composition 
of exponentials of the type exp(ai(t)Vi), where ai for i = 1,2, ... ,ci are scalar functions. In 
this manner it becomes possible to increase the order of the approximation without increasing 
the number of exponentials to evaluate and multiply together. We study order conditions and 
implementation details and conclude the paper with some numerical experiments. 

1 Introduction 

Although numerical methods for the integration of ordinary differential equations (ODEs) based on 
the use of the matrix exponential have long history, the subject has acquired new relevance recently 
with two developments. The first, which is irrelevant to the theme of this paper, is the introduction 
of Krylov subspace techniques and their application to large stiff systems of differential equations 
(Hochbruck, Lubich & Selhofer 1998). The other development is motivated by the philosophy 
of geometric integration and its purpose is to recover under discretization important qualitative 
and geometric features of the underlying dynamical system. Examples of such methods can be 
found inter alia in (Casas 1996, Crouch & Grossman 1993). An important technique in geometric 
integration is the use of Lie-group actions, which lend themselves to the design of very effective 
time-stepping methods for ODEs evolving on homogeneous manifolds. Such methods have been 
recently studied in (Munthe-Kaas 1997) and (Eng0 1998). Methods based on the use of the classical 
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Magnus and Fer expansions for integrating ODEs on Lie-groups can be brought into this formahsm 
(Iserles & N0rsett 1997, Iserlcs & N0rsctt 1999, Zanna 1997). AU such methods require a repeated 
evaluation of a matrix exponential, often of large matrices. Inasmuch as typically one can expect the 
replacement of the exact exponential by a suitable approximant (a rational function, say, a Krylov 
subspace approximant or a Schur factorization), the context of Lie-group methods imposes a crucial 
extra requirement. The approximant in question, applied to an arbitrary element of the Lie algebra 
0, must produce an outcome in the Lie group G, otherwise the whole purpose of the calculation, 
dicrctizing within G, will be null and void. This can be done is some, but by no means, all Lie 
algebras of interest and we refer the reader to (Celledoni & Iserles 1998) for a more substantive 
discussion of this issue. 

Let G be a finite-dimensional Lie group. For all practical purposes, we may assume that G is a 
subgroup of the general linear group GL(n), the set of all nonsingular nx n matrices. We denote by 
g the Lie algebra corresponding to G, observing that it is a subalgebra of 0l(n), the Lie algebra of 
all n X n matrices. Our concern in this paper is with differential equations that evolve on a manifold 
M subject to the action of G. For simplicity we can assume that M coincides with G and the action 
is of G on itself. The numerical solution of such differential equations can be obtained considering 
the pull-back on g, by means of the exponential map, of the vectorfield defining the equation. 
We can compute the corresponding flow by a Lie-algebra discretization method and recover the 
approximation of the original problem via exponentiation. Given an integration method of order 
p, we consider order-p approximants F{tB) for cxp(ii?), where B ^ Q and t > 0. We require that 
F{tB) G G, whence it is easy to prove that important qualitative features of the original equation 
and the order of the discretization are retained. In (Celledoni & Iserles 1998) we have introduced 
low-rank splitting methods for the construction of the approximant as the first attempt to provide 
a comprehensive treatment of this issue. 

Although the constraint F(tB) £ G represents remarkable advantage in many applications, such 
as problems in which the conservation of invariants is at issue in numerical modelling (volume 
conservation in meteorology, invariance under rotations in the theory of mechanical systems and in 
robotics), it should not be interpreted as the sole purpose of our analysis. Our methods are relevant 
also for the approximation of exp(<i3) in the more general setting B 6 Q\-{n). Suppose in fact that 
B G Q\-{n) and we want to approximate exp(i_B). It is always possible to write i? as a sum of a 
matrix Bg £ s[(n) (the special linear algebra oi n x n matrices with zero trace) and a diagonal 
matrix Bd whose nonzero entries are equal to 5 = tr {B)/n. Then Bg = B — Bd and [Bg, Bj] = 
so that exp(ii?) = cxp(ti?s) exp(M). This fact is a particular case of what is known in Lie theory 
as the Levi decomposition (Humphreys 1972, Varadarajan 1984). Using this decomposition of the 
matrix B, if necessary in tandem with some scaling and squaring technique, the approximation of 
exp{tB) can be always reduced to the approximation of exp{tBs) with Bg G s[(n). As long as we 
can assure that our approximation of exp^tBg) resides in SL(n), the outcome is an approximant 
F{tB) of exp(ti?) that shares with the exact exponential the feature that AetF{tB) = exp(tri?). 

It is possible to prove that, given a splitting B ~ X^iLi function 

known as the generalized Strang splitting, approximates exp(ii3) to order 2. As long as Bi, B2, ■ ■ ■ , Bk G 
g, it follows at once from the definition of a Lie group that the approximant resides in G. Moreover, 
2fc — 1 is the least number of exponentials that render such a splitting into a second-order approx- 
imant (Celledoni & Iserles 1998). The Strang splitting is time reversible, hence it follows readily 
from classical theory that the order can be raised from 2 to 4 by composing three Strang splittings 
with different time steps (Yoshida 1990). In that case we need to evaluate 3fc exponentials and 
multiply 6k matrices. In the case of low-rank splittings which have been considered by Celledoni & 
Iserles (1998) this results in the following count of flops: 4n^ for order 2, 12n^ for order 4. 

In this paper we present composition methods in which the number of exponentials k equals the 
dimension d of the Lie algebra. Our construction allows us to increase the order of the approxi- 
mation without increasing the number of exponentials to evaluate and multiply together. Letting 
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{Vi, . . . , Vd} be a basis of g, we write F{tB) as a composition of exponentials of the type cxp(ai{t)Vi), 
where each ai{t) for i = 1, . . . , d is a scalar function. In general d = 0(n^) , however, with an appro- 
priate choice of the basis elements, the computation of each exponential exp{ai(t)Vi) requires 0{n) 
flops, while the formation of their product adds just 2n^ + 0(n?) flops. The challenging part of the 
computation is the construction of the functions ai(t), . . . ,ad{t), and the cost of their calculation 
depends on the desired order of the approximation. Naive complexity analysis might have indicated 
that the total cost is growing exponentially in d as the order increases. Yet, the cost remains rela- 
tively modest for small orders and the method lends itself very well to the exploitation of sparsity 
in the matrix B. In the sequel we show how this approach can be turned into an efficient numerical 
method and we obtain algorithms of order up to 4 with a cost of 0(n^) for dense matrices. 

Our approach can be interpreted as representing the solution using canonical coordinates of the 
second kind, an approach that has been pioneered by Owren & Marthinsen (1998) in the context 
of general Lie-group methods. Having said this, the more restrictive framework of exponential 
approximants possesses a very great deal of special structure. This can be exploited so as to 
produce efficient and competitive algorithms that approximate exp(ti?), i? G g, in the Lie group G. 

2 The technique of coordinates of second kind for the ap- 
proximation of the exponential matrix 

Let G be a Lie group and g its corresponding d-dimensional Lie algebra. We choose a basis 
{Vi, . . . ,Vd} of g, whence every element F S G sufficiently close to the identity can be represented 
in a unique fashion as 

Y = exp(7iyi) exp(72T/2) • • ■ exp(7rfVd), 

where exp : g — > G is the exponential map. This representation is known as representation in 
canonical coordinates of the second kind (Varadarajan 1984). This representation is global in the 
case of solvable Lie algebras. We restrict ourselves to the case C gl{n), G C GL{7i), when exp is 
the usual matrix exponential. 

Given _B S 0, we can represent it in a unique fashion as 

d 
i=l 

It is possible then to write cxp{tB) in the form 

C/(t) = exp(iB) = exp(gi(t)l^i) exp(.g2(t)l/2) • • ■ cM9dit)Vd). 

Letting g = [gi, . . . , gdf^ , (3 = [(3i, . . . , Pd]"^ , it can be proved that the vector function g obeys a 
differential equation of the form 

^ = f{f3,9), 9(0) = 0, 

where / is a suitable function of (3 and g, for sufficiently small t (Wei & Norman 1963). Given 
a solvable Lie algebra Wei & Norman (1963) prove results on the global representation of U. 
However, an explicit form of / is known only for very simple examples of low-dimensional Lie 
algebras. 

In this paper we seek polynomials ai gi, . . . , ad ~ gd of a suitable degree so that 
exp{tB) « exp{ai{t)Vi) exp(a3(i)V2) • • • exp{ad{t)Vd) . 

Differentiation yields 

j2m = J29'.it)U^'''^'y^ tl (2-i) 
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Evaluating this expression at the origin gives the first-order condition 

5:(0)=ft, z=l,2,...,d. (2.2) 

Further differentiations of (2.f ) lead to higher-order conditions. Let us define the functions 

Pi(g) = exp(adgiVi) o ■ ■ • o cxp(adg,_^VA„J(yi), i^l,2,...,d, (2.3) 

where the adjoint operator ad^; : g ^ g is defined as adxiy) = [a;, y] for any x,y € [x, y] = xy — yx 
being the matrix commutator. Note that Pi{g{jS)) = Vi, i = 1, 2, . . . , d. Moreover, the right-hand 
side of (2.1) can be written in the simplified form The function 

d 
i=l 

Since the derivatives of the left hand side of (2.1) vanish, the conditions for order p > I, can be 
obtained by solving the equations 



0, r==l,2,...,p-l, p>l, (2.4) 



where 



In particular. 



d 



i=l k=l 



irto,=i:(,rp.+4p.,. 

1=1 ^ 



dt2 ^ V" ' dt ' ^Mt2 

i — 1 ^ 

i=l ^ 

Solving (2.4) for r ~ 1 results in the values of gf (0) for i = 1, 2, . . . , d that allow us to construct 
an order-2 approximant. Substituting such values in (2.4) for r = 2 yields g'/'{0) for z = 1, 2, . . . , c? 
and consequently an approximant of order 3. Similar procedure can be used to construct recursively 
approximants of arbitrarily high order. 

The main part of the computation is the evaluation of the fc-th derivative of Pi{g) at t = 0. 
Expanding the exponentials in (2.3) we obtain 

i-l 

P^{9) + ad,,y, + iad^^^^ + iad^^^^ -t- . . .)(F,) 

k=l 

and, after further algebra, 

{i-l i-l fe-1 

I + Y^ adg^v^ + E E ^dgjV^adg^y, 
k=l k=2 1=1 

i-l i-l fe-1 (-1 

+ m ^^Lvfc + E ads, v-^ adg,v,adg,y, 

fe=l k=3 1=2 3 = 1 

i-l fe-1 i-l "I 

+ 5 E E i^'^a.yMlv, + ad^,^,adg,yj + i ^^Iv^ + • ■ • (^0 • 

k=2 1=1 k=l ) 
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Similarly to (Owrcn & Marthinscn 1997), we write Pi{g) in the form 

oo i—1 2 — 1 i—1 



OO i—1 2 — 1 i—1 ^ 

Pi (g) = / + ^ ^ ^ ■ • • Tjfii • ■ • 9jMv,, o ■ • • o ady^^ {V,). 



r = ljl = lj2=jl jr=jr-l 

Here j = (ji, . . . jr) is a multi-index of integer elements with 1 < jr < * ^ 1 and j! := (ji!(j'2! • • ■ Qi-i^- 
where qk is the number of occurrences of k in {ji,j2 ■ ■ - jr)- 

A general expression for the fc-th derivative of Pi is given as follows: since gi{0) = 0, we may let 
fi{t) = gi{t)/t, i = 1,2, . . . ,d. We can then rewrite Pi in the form 

oc ? — 1 ? — 1 i—1 

= / + E E E ■ • • E 71 • • • /^-ady^i ° • • • ° adv,. (K). 

'" = 1 jl = lj2=jl jr=jr-l 

By following the construction in(Owrcn & Marthinscn 1997) we obtain 



d''P^ 



t=0 



E E 

r=l Si + ...+5r = k 



k\ 



E 



<jV,<'-i ■ 



(2.6) 



t=0 



ady^^ o • ■ • o adv^^ {Vi). 



Substituting (2.6) in (2.4) and (2.5) we obtain the conditions for arbitrary order p. In particular 
we obtain the following formulae for the derivatives of Pi{g) aX t = 0, 



dp 
At 

d^P^ 



dt2 
^?P^ 



di3 



t=0 



k=l 

i-1 /fe-1 \ 

E E 2adv,ady,(K.)5U0)5;(0) + ad^^ (F.)[5l.(0)]2 + ady, (K^glf (0) , 
k=i \i=i / 

i-l k-1 /-I 

6 E E E ^dv. ady, ady, (KOg; (0)gK0)5fe (0) 
fc=3 ;=2 j=i 
i-i fc-i 

+ 3 E E (adv^ad^.^ (T^.)<?KO)[<?UO)]' + ad^^-^ady, m[9m?9km) 

k=2 1 = 1 
i-l 



+ Ead^,(1^0[5fe(0)]' 

fc=i 

i-l fc-l 

+ 3 E E adv.ady, (V,) (.gj:' (0)51(0) + g'M9im 

k=2 1=1 
i-l 

+ 3Ead'yjy,)5fe(0)5fe(0) 

k=l 
i-l 

+ Y,aAvM)9k{Q)- 



k=l 

Substitution readily produces order conditions. Specifically, 

d d i—1 



Y.9iim = -E5K0)Eady.(K)5l.(0), 



(2.7) 



i=l 



k=l 
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are conditions for order 2, while 



k=l 
k-1 



1=1 



fe=l 

+ ad2y,(K)(gU0))' + adv,(F,)3fc 
are the order-3 conditions. Finally, conditions for order 4 are 

d d ( i—1 



(2., 



i=l 1=1 I 

+ 3gr(0) 



fc=i 
i-i /fc-i 



^ ^2adv,ady,(F,K(0)ft'(0) 

ad'y.(^^«)(5fe(0))' + adv,(F.K:(0))] 

j-i fc-i ;-i 

.9K0) 6 E E E adv^,ad^,ady, (K.)<?;- (0)5^(0)5^0) 



fc=3 ;=2 j=i 



i-l fe-1 



+ 3E E (adv,ad^jy05;(0)(5U0))' + ad2v-,ady,(V;)(5K0))25U0)) (2.9) 
fe=2 ;=i 

k=l 

i-l k-1 

+ 3EEadv,advjy.) (5l^(0)5;(0) +5U0)5r(0)) 
fe=2 ;=i 

+ 3Ead'y.(^05fe (0)5^0) 

k=l 

+ Y.advM)9kiO) ]■ 

k=l J J 

In Figure 1 we have plotted along the y axis the 2-norm of the error of the approximation of 
ex-p{tB) with the second-kind coordinates (SKC) methods of order ranging from 1 to 4. The values 
of the error are plotted against time, (along the x-axis), to logarithmic scale for matrices of s[(5). 
The methods have been implemented using the standard basis defined in section 3. 

The computation of 5"(0), 5"'(0) and 5^^(0) is obtained directly implementing the formulas 
(2.7), (2. 8), (2. 9) respectively. This implementation does not depend on the choice of the particular 
basis of 5l(5), but the number of commutators that must be computed with this approach is 0{dP) 
for p = 2,3, 4. Even if we assume that the ViS are very sparse matrices and that the cost of 
computing each commutator is 0(1) operations, the total cost exceeds 0{n^P) flops for p ~ 2,3,4 
where n is the dimension of the matrix. Such expense is not acceptable for a competitive method 
of approximation of exp(ti3). Fortunately, it can be decreased a very great deal by an appropriate 
choice of the basis {Vi, V2, . . . , V^}. This is the theme of the next section. 
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SKC for iraee-free matrices n=5, order one to four 




Figure 1: Error in the approximation of the exponential with WN technique. 



3 Choosing a basis 

The choice of the right basis and sparse representation of commutators are critical to the implemen- 
tation of the SKC methods. Recalling the order conditions (2.7)-(2.9), our aim is to choose a basis 
so that terms of the form a.dvi^ ady^^ • • • ady^^ Vj can be represented in the most economical man- 
ner. We recall that, given the basis {Vi, V2, . . . , Vd} of a d-dimensional Lie algebra g, the structure 
constants arc the numbers k,l,i = 1, 2, . . . , d, such that 



{Vk,Vi]=J2cf,^,V, 



i=i 



(Humphreys 1972). Let 



k=l 



Then an order- 1 condition is always 

.9fc(0)=A, 



fc = L2, 



(3.1) 



To obtain the order-2 condition we substitute (3.1) in (2.7) and express commutators in terms of 
structure constants, 

E aum = - E A E [^j- ' = - E A E E sv^^ 

fe=l 1=1 j=l 1=1 j = l k=l 

d / d l-l \ 

= -E EE/^^sV.V- 
k=i j=i / 

Since cjj = — cf we thus deduce that 



d l-l 

1=1 j=i 



(3.2) 
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Likewise, substituting in (2.8), 

d d ( i-l i-l I l-l 

E 9k im = - E i 25r(o) Et^' ^^f^^ + a E hEt^^' iVi^^MPj 

k=l i=l [ 1=1 1=1 \ 3 = 1 



Note that 



d d 



[V,, [Vu V,]] = E ^A^i^^A = E E <»s>^fc- 



s=l k=l s=l 



Therefore 



E - -2 E 5f (0) E E -''-ilP^Y^Y.ilil ^lA^^m^^ 

fc^l j^l k^l s^l 

d i—1 d d d i—1 d 

-Y.p^T.T.Y. <As0fy^ - E E E <^yamk 



i=l 1=1 k=l s=l 1=1 1 = 1 k=l 



and we deduce that 



d i-l d i-l l-l d 

9k {0) -EE <i\^9-mi + P^9lm + 2E E E E <i4sPM 

i^l 1^1 i^l 1^1 j^l s^l 
d d i—1 d 

-EEEE<^^t/3»A', fc = i,2,...,d. 

1=1 i=l 1=1 s=l 

Bearing in mind that for order p we require 

p 



W = E45i'''(0r, A: = l,2,...,d, 



^ r! 

r=l 

we observe that the sheer volume of calculations required for the evaluation of the functions 
«!, a2, . . . , is prohibitive for, say, order 3, unless most of the structure constants vanish. For- 
tunately, bases of finite-dimensional Lie algebras which are 'sparse' (in the sense that a very high 
proportion of structure constants vanish) are known. They are associated with root space decom- 
positions of Lie algebras (Humphreys 1972) and, in the case of semisimple algebras, are known 
as Chevalley bases (Carter, Segal & Macdonald 1995). Wishing to avoid too much Lie-algebraic 
terminology in a numerical analysis paper, we reserve our exposition to just three examples which 
are the most important in a range applications. 

The orthogonal group Let Q = 50 (n), the Lie algebra of ti x n skew-symmetric matrices. It cor- 
responds to two important Lie groups: the orthogonal group 0(n) of n x n orthogonal matrices and 
its subgroup, the special orthogonal group SO(n) of matrices with unit determinant. Its dimension 
is c? = ^n{n — 1). We let 

Fij ~ EiBj^ ~ EjEi^ , i = l,2, ...,n, J = i + 1, i + 2, . . . , n, 

where is the z-th canonical vector of M". In other words, Fi^j is a matrix whose (z, j)-th element 
is 1, the (j, i)-th element equals —1 and zero otherwise. We can trivially expand each B G SO(n) as 
B ~ '}2^=iJ2^j=i+i^i,j-^i,j- U{t) :— cx.]}{tFi j) is simply an Euler rotation in the (i,j) plane: it is 
identity matrix, except that 





U^.3 ' 




cos{tFij) 


siii{tF,j) 


. U3,r 


Uj,3 . 




-sin(tFij) 


cos{tFij) 
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Noting that 



— Fid, 
Fi,k, 

O, 



i 7^ ^ j = k, 
i^k, j = 1, 
i = k, j ^l, 
otherwise, 



the order conditions arc simpHfied as follows, 
p>l: <,(0) = 6„ 

P>2: 5-:,(0) 



r—i+l 

1 < i < j < rt. 



and similarly for higher-order terms. Thus, the cost of computing the coefhcients for the second- 
order method is just ^{n — 2){n — l)n « ^n^ flops. In comparison, a naive computation of (3.2), 
without exploiting sparsity of structure constants, requires ^(n^ — n — 2)(n — 1) 



2^2 



1^6 



rr flops. 



A more classical composition method for B £ 50{n) is the Strang splitting which we can write 
in the form 



atfel, 2^1,2/2 , , , tb„ 



^Fn-l,n^tb„ 



-2,>>/2 . . . pt6l,2-Fl,2/2 



(Celledoni & Iserles 1998). It gives a second-order approximant to exp(t_B) whose calculation re- 
quires w An^ flops, in comparison with w 3n'^ for the second-order CSK method. 

We note that, in the specific case of SO(n), diagonal Pade approximants to the exponential 
provide an alternative to our method, since they map the algebra to 0{n). Having said this, for 
dense matrices B the cost of evaluating the second-order approximant (/ — ^tB)~^{I + ■T^tB) with, 
say, LU factorization is 0(n"^), comparative with our method. 

The special linear group Let = sl{n), the set oi n x n matrices with zero trace, whence 
d = — 1. We split the algebra in the first instance into diagonal and off-diagonal parts: in the 
terminology of Lie algebras, the subspace spanned by the diagonal elements is a Cartan subalgebra 
(Carter et al. 1995) or maximal toral algebra (Humphreys 1972) of s[(n). Specifically, our basis is 



«,J = 1,2,. 



7^.?}U{A : z = 1,2,. 



where 



Di = BiBi B.i-^iB.^^i, 



z = l,2,. 



1. 



The exponentials of Eij and Di are trivial. 



I + tEi 



We order the elements by taking first Ei^j , i j, in lexicographic order, followed by Di , D2, • ■ ■ , F)„ 



The commutator table is 



r E, 



[Fi,j,Er^s\ 



t.s 1 
~Fr,j, 



i = s, 3 7^ ^ 
i = s < j = r, 
i ~ s > j = r, 
otherwise. 



hj,r,s = 1,2,. 
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[Eij , Dr] — 



( p 


i^r, j + 


1, 




■i'^r, j = r + 


1, 


— 2£'r,r+l, 


i = r, j = r + 


1, 


Er+l,j, 


i = r + I, j =/= 




Ei_r 1 


j = 


r, 


2£'r+l,r, 


i = r+l, j = 


r, 


o, 


otherwise, 





ij = l,...,n, i^j, 
r = 1, . . . , n — 1, 



[D,,D,]^0, = 1,2,..., 71,-1. 

In general, for a d-dimensional Lie algebra there are {d — structure constants. In the case of 
s[(n) this means that up to « nP structure constants may be nonzero. Yet, using the above basis 
results in just 2(n — l)n^ + 4(n — 2)(n — 1) + ^{n^ — l)n « ^n^ nonzero structure constants and 
substantial saving in the implementation of the SKC technique. 
Letting 



(kj) k 
[Eij,Dr] = ^ C(^ij^^rEk,l 

(note that [Dr,Eij] = —[Eij^D,.] and [Di,Dj] = O) we thus have 



(fe,i) 



+1, 


k = i, I ~ 


s, r^j, i, 


-1, 


k = r, I = 


j, r j, s = i, 


0, 


otherwise. 




+1, 


i = s <j - 


= r, fc = s, s + 1, . . . , 7' — 1 


-1, 


i^s> j - 


= r, fc = r, r + l,...,s — 1 


0, 


otherwise. 




+1, 


k = i = r - 


f 1, I = j, j =/= r OT k ~ i, 


-1, 


k = i = r, 


^ = .h + 1 or k~i, 


+2, 


k = i = r - 


f 1, l=J=r, 


-2, 


k = i = r, 


/ = j = r + 1, 


0, 


otherwise. 





„(k,l) ^^k ^ Q_ 



Letting 



B = Y,Pk,iEk,i +Y.^kDk, 

k^l k 



and ordering the pairs (fc, I), k I, in lexicographic order, we thus have 

9'kM = Pk,i. 
9'kiO) = lk, 



9k.i{0) 



{i.j)y{r,s) {'i,j),r 
k— 1 n 

E Pk^iPiJ - E Pk,il3i,l + Pk,l{lk-l + 7i - 7fc - ll-l), 



i=fc+l 
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where 7o = 7n = 0. 



The Lorenz group This is the 6-dimensional group SO (3, 1) of 4 x 4 matrices A such that AJA^ = 
J, where J = diag(l, 1, 1, —1) (Carter at al. 1995). It has important apphcations in special relativity 
theory, he corresponding Lorenz algebra SO (3, 1) consists of all matrices B such that BJ+JB"^ = O. 
It is easy to verify that each clement of S0(3, 1) can be written in the form 



B 





-b2 
bs 



hi 62 &3 

&4 &5 

-64 66 

&5 &6 



61,62, ... ,65 e 



Choosing the basis 





0100" 




0010" 




"0 


000" 




"0001" 




"0000" 




"0000" 






-10 




000 







1 









1 









[ 





7 


-1000 







-10 














1 


> 






















10 




10 




10 





we obtain the commutator table 
















[^1,1^2] --"t^3, 




V2. 












= 0, 


[^2,1^1] =^3, 






[V2M 








[V2,Vq] 




[^3,1/1] =-^2, 


[^^3,^^"2] = 






= 0, 




= -v&. 




= ^5, 


[Vi,V{\=V5, 


[Vi, V2] = 






= 0, 


[^4,^5] 






^V2, 


[V5,V{\ = ~Vi, 


[^5,^2] = 


0, 


[^5,^3] 


= ^6, 






[^5,^6] 


= V^3, 




[^6,^2] = 




["^^6,^3] 


= ~V5, 




= -V2, 




= -F3 



Thus, out of 180 structure constants, just 24 are nonzero - and they all equal ±1. After brief 
claculation, we drive for example the polynomials that yield an order-2 CSK approximant, 



ai{t) 
a2{t) 

ai{t) 

aQ{t) 



Pit- 
Pst- 



i(/3i/33 
i(/3i/32 



P4P5)t\ 

Pipe)t^, 



Pit - Ufiifi^ + P2P&)t^ 



P5t- 

Pet- 



\{PiPa 

l(/33/35 



P3P^)t\ 

P2Pi)t', 



where B = Y,\^^PkVk. 



4 Time symmetry 

An approximant F{tB) w cxp(t_B) is said to be time symmetric if F{tB)F{--tB) ^ I, t > 0. Time 
symmetric approximants are important for a number of reasons, not least being that they lend 
themselves to the Yosida technique, which allows their order to be increased (Yoshida 1990). The 
techniques of the last section are not time symmetric. Here we describe their modification, which 
results in a time-symmetric approximant. 

Me mention in passing that it is possible to envisage two distinct techniques to obtain high-order 
algorithms based on canonical coordinates of the second kind. The first, implicit in the work of the 
previous section, consists of evaluating the numbers (7^'^(0) for Z = 1, 2, . . . ,p, where p is the order 
of the method. The alternative, the subject matter of the present section, consists in combining 
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a second-order or a fourth-order approximant across a number of steps to obtain a higher-order 
method. 

Given the sphtting 

1=1 

it is weU known that the Strang splitting The approximation 

FitB) - e*^i/2 . . . gtc._i/2gtc.gtc,_i/2 . . . gtCi/2 (4^^) 

is of order 2 and time symmetric. Note that, as a consequence of time symmetry, for sufficiently 
small i > we can represent F{tB) = e^^*) where the matrix function !F{t) is odd. It is precisely 
this feature that allows the application of the Yosida technique. 

The clear reason for (4.1) being time symmetric is that it is palindromic in the alphabet 
{Ci, C27 • ■ ■ , Cs\- This provides a clue how to modify techniques based on canonical coordinates 
of the second kind so as to render them time symmetric. Given a basis {Vi, V2, . . . , Vd] of the Lie 
algebra 0, we approximate e*^ by the product 

exp[ai(t)Vi] • • •exp[ad_i(i)Vd-i] exp[ad(01^d] exp[ad_i(t)Vd_i] • • • exp[ai(t)Vi], (4.2) 

where ai, 012, . . . , are odd polynomials. 

Taking ai = \(3it, I = 1, 2, . . . , d — 1 and ad = Pdt yeilds the second-order Strang splitting. In 
the sequel we seek higher-order methods of this kind. 

Using the Baker-Campbell-Hausdorff (BCH) formula it is possible to express the product of 
exponentials at the right hand side of (4.2) as a single exponential (Varadarajan 1984, p. 141). 
Due to the symmetric arrangement of the exponentials in (4.2), the BCH formula is an expansion 
in odd powers of t. If this expansion converges, which is always the case for sufficently small t, it 
makes sense to write the equation 

d—l 00 

tB = 2Y,a,{t)V, + ad{t)Vd + ^ Q'^(a). (4.3) 

i=l k=l 

Here we denote by Q'^'^{oi.) the terms of order C'(t^'^+^) in the BCH formula apphed to (4.3). 
Moreover, we let af'it) be the polynomial obtained by truncating the expansion of ai{t) after the 
first k terms, and we denote the remainder by rf^{t). In other words, 

a,{t)=af{t)+rf{t), rf (i) = , z = l,2,...,d. 

From (4.3) we deduce 

2Y,af''\t)V, + aT\t)Vd = tB - J^Q'^i^) + ^l^'"^') • 

i—1 r=l 

Noting that 

Q^^ia) = g2'-(a2("-i) + r^C^-D) = Q2r(^2(fc-i)) ^ o^^2k+r^^ ^ 

we obtain 

2j2c^f'\m + aT\t)Vd = tB - Q2.(^2(fc-i)) ^ 0(^2fe+i) (4 4) 

i—1 r—1 

Dropping the C'(t^'"'+^) terms in (4.4), it is possible to compute a^*' from a^C^^^'. This gives a 
procedure to derive a sequence of successively increaing-order approximants of exp(ti?). It is easy 
to see that the approximants 

F^-f^itB) = exp(af (t)l/i) • ■ • exp(af (t)^^) ■ • • exp(af (t)l/i). 
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of cxp{tB) are such that F^'' (tB)F^'' {—tB) = /, hence time symmetry, the reason being the sym- 
metric arrangements of the exponentials in F^^{tB) and the odd-power expansion of the functions 

The BCH and symmetric BCH formulae for fc-terms have an exceedingly complicated expansion, 
which can be obtained recursively. In what follows we will make use just of the term (3^(a), 
demonstrating how it is possible to compute it explicitely for particular choices of the basis. 

In the remainder of this section we consider the implementation of time-symmetric CSK methods. 
We split B as before and commence by considering the Strang splitting (4.1) except that, to simplify 
notation, we arrange the terms in reverse ordering. 



Lemma 1 The term of the BCH formula applied to (4-5) is 



= ^ Y.[Ci + ■ ■ ■ + + iC,, [Ci + ■ • ■ + G] 



(4.5) 



(4.6) 



1=2 



Proof See the appendix. 



Let us next consider the case = 50(n), choosing the same sparse basis as in Section 2. Therefore, 
according to (4.6), we have 



i=l j=i+l 



(4.7) 



i— 1 j=i+l 



We compute separately each part of this sum. Exploiting the commutator table of our basis we 
have 



[^1, 2^^1,2 + . . . + bij^iFij^i,bijFij] 

/ j-i i-i i-i 

bi,j[- H bj^sFs,] - yi brjFrs + yi bt,iFt^ 



3 ' 



and noting that bi^g = —bs^i, we deduce that 

[^1,2-^1,2 + . . • + feiJ-l-PtJ-l, 6ijfi.j] = b^j 

Commuting the right-hand side with Fij gives 



i-l 



t=l 



r=l 



Let 6i , 62 , ■ • ■ , be the columns of B and denote 

s-l 

bi = Hfcfe.zefc, 
fc=i 



f 3-1 i-l 

H bt.iFtJ — H brjFr^, 
1 t=l r=l 



J 



i-l i-l 

yi bt,iFt^i + H brjFrJ. 

t=l '-=1 



fc = 1, 2, . . . , n. 



(4.8) 



(4.9) 
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Then (4.8) yields 

[foi,2i^i,2 + . . . + 6,,j.F,j] = h,j{biej - Bjbf) - bi^j{b)ej - e,bf), 

while (4.9) gives 

[Fijj [^1,2-^1,2 + ■ • • + bij^iFij^i, bijFij]] 
= b,,,(Mef - e,bf ) - 6,,,(b}eJ - e,bf ). 

Multiplying the latter by bij and summing in i and j we can evaluate (4.7) in operations. Note 
that we count separately multiplications and additions, for example, we assume that the cost of 
Euclidean inner product of two vectors of length n is 2n operations. 

We now assemble together our results to calculate (4.7). We proceed by splitting the sum 
bi,2Fi^2 + ■ ■ ■ + bij-iFij-i in three parts, whereby 

[6i,2^\,2 + • • • + b^ej - e,6f - {b]ej - e.bf )] 

= E EK^i^T _ e,6f ), b^ej - e,lf - (b^^ef - e,bf)] 

1=1 k=i 

m i—1 

+ E E[(^?^^ - - e,bf - {b]ej - e.bf )]. 

Finally, 

[6i,2i^i,2 + ■ • • + 6^eJ - ejb(^ - {b]ej - e,bf)] 

i-l 

biMeJ - eM^) - b^i^,, - bi^,{b\ej - e,bf ) + bf b}^],, 

1=1 

5] 5,,(b?+ieJ - e,bj+i^) - bj+i^b^i^,, + b,,ib)e[ e,bf ) + b^+^^bjf), 

n 

+ J2 lfb)Fu~lfb^Fi,. 

1=3 + 1 

We analyse the computational costs of the previous formula, summing over i and j and showing 
that (4.7) can be computed in O(n^) operations. Note that, since X^tZ+i bi,i^i = ~ have 

b^bf - b}bf + ^ 6u(b}eT _ e,bf ) - (bjbf - b^bf ) = -2(b}bf - b^b}^). 

i=4+l 

It is more convenient to write the previous expression in the form 

i—1 
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n / n 



n — 1 n 

j 



J2 E 2h,(b^m~bri'^{bi-bDh 

i=l j=i+l 



The first part of this sum is computed in about operations and the second part, exploiting the 
equality 

n-l n n-1 i+2 /k+l \ 

i=l j=i+l 1=1 k=n—l \l=n / 

can also be computed in |n'^ operations. 

Adding terms of the type cxFij and leads to 





1=1 


E ^'^KFi, - 

l=i+l 


n 

E ^'^b^^). 


and 


Eferfe}^.- 

1=1 




- E ^^-^ 
i=j+i 



where the matrix Li is the lower triangular part of 6i.2-F'i,2 + . . . + 6i_i_„i^i_i_„ and we denote by 
L|, s = the matrix Li with zeros along its s-th row. Summing up with respect to i and j, we 
obtain 

n— 1 n n—1 / n \ 

n — 1 n n / \ 

E E f'^^jFibjej^ = ^ E^^.j'^* ) 

i=l j=i+l j=2 \i=l 



T 



where we have used the notation C; := Lib] for i = — 1. The cost of computing the first sum 

I' 



is In'^, while the cost of computing the second is 2n^ operations. 



Finally the terms 

n—1 n i n—1 / n \ / ^ \ 

E E E^^.^M(b^ej - e.-bf) E ''M E ^.^'-J - E 

i=l ]=i+l 1 = 1 1=1 \j=l+2 / \j=l+2 J 

with Cj,i = bi.jbi^i, 

n—1 n i n—1 / n — 1 \ / n — 1 \ 

i=l j=i+l 1=1 1=1 \i=l+l ) \i=;+l / 

with di^i = I]"=j+i bijbij, and 

n—1 n j— 1 

E E E b,,k,{l^'ej-e,br^) 

i^l j^z+1 l^i-\-l 
n — 1 n—1 

= E E 6M(br(e^-b.r-(br^-bOb: 



i=i i=i+i 



i+l J 
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can be computed in about |n'^, |n'^ and ^n^ operations respectively. Collecting the contributions 
of all the terms in the sum we obtain a total count of 7^n^ operations. 

At the present time it is not clear that this method of computation of in the SO(n) case is 
optimal form the point of view of complexity theory. We did not try any other ordering of the basis 
elements and it is not at all certain that different orderings could give better constants in front of 
the term n^. 

Given that the construction of the (second-order) Strang splitting carries a cost of An^ oper- 
ations, the total flop count for constructing a symmetric fourth order SKC approximation of an 
exponential in SO(n) by our algorithm is 11 ^n'^. This is marginally better than obtaining an order- 
4 approximation by the Yosida technique from three Strang splittings which, as pointed out in 
(Celledoni & Iserles 1998), requires 12?i^ flops. 



5 Sparse matrices 

In a naive formulation, the method of canonical coordinates of the second kind is considerably too 
expensive for practical computation. This, however, can be alleviated by the use of a sufficiently 
'sparse' basis of the underlying Lie algebra Q. As explained in Section 2, choosing a basis so that an 
overwhelming majority of structure constants vanish renders the algorithm strikingly more effective. 
It is important to emphasize that this has nothing to do with the structure of the matrix _B G g, 
which need not be sparse. Yet, in most practical computations (in particular when n is large) 
one can expect B to be sparse and structured. Good algorithms should be able to exploit this 
phenomenon. 

In the case of SKC methods wc identify two mechanisms that allow us to exploit sparsity. 
Although this aspect of our methods is still a matter for active investigation, the interim results are 
substantive enough to warrant publication. For simplicity, we describe the first mechanism just in 
the case of a tridiagonal B G SO (n) . hence 



B 



Pi 
-Pi 










-Pn 



" 







n-1 

— ^ PkFkM+1 




fe=l 


Pn-1 




1 





where the matrices Fk.i = SksJ — eiej have been introduced in Section 2. Since 



Pk, 
0, 



l^k + l, 
otherwise, 



it is easy to substitute in the general formulae for the ordcr-2 method: 

n' (0) = I J=« + l, q"(0) = I 2, l<i<j<n 

^^'^^ ' I 0, otherwise, ^'■'^ ' \ 0, otherwise, ~ J — ^ 



Arranging the elements of the basis in lexicographic order, we thus obtain the second-order approx- 
imant 



In other words, the cost of the approximation is just 0{n) flops. 
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Similar situation pertains to 



B 



71 Vi 
Ml 72 ' ■ ■ 

■■• ■■• 











In 





Vn-l 

In 



e 5\{n). 



Choosing the same basis and terminology as in Section 2 we can readily ascertain that 

5fc,fe~2(0) =Mfc-2Mfc-l fc = 3,4, ...,71 

5fe,fe-i(0) = -(7fc-2 - 27fe_i +7fc)^fe_i, fc = 3,4, . . .,n, 
fffc,fe+i(0) = (7fc-i - 27fc +7fe+i)77fe, = 2,3, . . . ,n - 1, 

fffc,fc+2(0) -rikVk+i, fc 1, 2, . . . , n - 2, 
<,(0)=0, |fc-Z|>3 

and <?Jj'(0) = —rjkfik, fc = 1, 2, . . . , n — 1. Thus, a second-order approximant to a tridiagonal B G sl(7i) 
is itself quindiagonal and its computation requires just 0{n) flops. 

Higher-order approximants and matrices with greater bandwidth lend themselves to similar 
treatment, although the savings are less striking. In a sense, the situation is parallel to that of 
approximating exp{tB) by a rational approximant, when savings accrue from sparse matrix-inversion 
methods, except that in our case the result is assured to belong to the right Lie group. 

Another observation which is highly pertinent to the approximation of exponentials of sparse 
matrices has been made in (Iserles 1999). Suppose that i? is a banded matrix of bandwidth s > 3. 
In general, F(t) = exp{tB) is a dense matrix. Yet, as is easy to illustrate by computer experiments, 
F{t) is very near to a banded matrix. Specifically, given e > 0, there exists r = r(t,e) > s such 
that all the elements of F{t) outside a band of width r are less than e in magnitude. Moreover, 
tight upper bounds on r can be derived with relative ease. The idea thus is to set to zero all the 
elements outside bandwidth r. The outcome is a banded approximant to the exponential. Moreover, 
with an appropriate choice of basis elements, this means that the functions ai are set to zero for 
elements that possess terms exclusively outside the band. Consequently corresponding exponentials 
equal identity and need not be included in the product. Thus, the cost scales with the size r of the 
bandwidth. Similar phenomenon has been already encountered in the context of 50{n) and si{n), 
when our choice of basis and order has implied a banded structure of the exponential. The present 
mechanism is different, even if the net outcome is similar. 



6 Numerical experiments 

Our numerical experiments are organized as follows. We fist consider a test on random matrices in 
SO (50), illustrating the performance of methods based on the use of second kind coordinates tech- 
niques for full and sparse matrices. The third and last example is the solution of a third-order ODE 
using Runge-Kutta/Munthe-Kaas (RK/MK) methods described in (Munthe-Kaas 1997). We use 
the Matlab toolbox Dif f Man for the integration of ODEs on manifolds, comparing the usual imple- 
mentation of RK/MK methods, whereby the the exponential is approximated to machine accuracy, 
with a version of the methods obtained using the time-symmetric fourth-order approximation from 
Section 4. 

All experiments have been performed in Matlab and we have computed the error while comparing 
the results with the built-in function expm which calculates the exponential to nearly machine 
accuracy. 
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SKC- symmetric order-4 and Strang splitting, skew- symmetric matrix n=50 




Figure 2: Error versus time in the so(50) (full case). 



We evaluated the the error computing \\e~^^ F{tB) — I\\f where F{tB) is the SKC approximation 
of exp(ti?) and || • ||f denotes the Frobcnius norm. The matrices have been generated randomly 
using the Matlab function rand and scaling the Frobcnius norm so that ||-B||f = 1- 

We approximate cxp(ti?) with a single step of the methods for different values of <, {t = 1/2^ 
and fc = 1, . . . , 5). 

In both the first two figures the norm of the error is plotted (along the y-axis) to a logarithmic 
scale with respect to t. Figure 2 reports the results of our first test, where we have considered a 
full matrix in so(50). In the plots the error norm is indicated with the symbols (SKC, time 
symmetric, order 4) and 'o' (Strang splitting, order 2). 

In the next example, illustrated in Figure 3, the same methods have been applied to a sparse 
matrix in 50 (50), with four non-zero diagonals (i.e., bandwidth 5). In both the examples the 
methods give the correct order. In the second case, however, the count of fiops is drastically 
reduced. We counted the number of fiops using the Matlab function flops. In the first case the 
cost for constructing amounts to 9.62n'^ while in the second we counted 0.95n'^ flops. As it is easy 
to understand, the described implementation of the methods allows to take advantage immediately 
of the sparsity structure of the matrix B, working directly on the nonzeros entries oi B a la Section 4. 

The last example is concerned with the use of the techniques described in this paper in substi- 
tuting the exponentials computed to machine accuracy in the integration methods of (Munthe-Kaas 
1997). The experiments have been performed using the Matlab toolbox Dif f Man . We use a RK/MK 
method of order four. 

The example is a problem whose solution is the soliton originating in the Korteweg-de Vries 
(KdV) equation. It is a third-order ODE obtained performing a symmetry reduction on the KdV 
equation. The resulting ODE can be written as a three-dimensional system. 



y' = 



1 
1 

-9y(2) 3 



with y(0) = [1,0,-1.5]"^ and t £ [0,5]. The solution of the ODE / = yi{t) can be easily derived 
explicitely and it is /(<) = asech(t/3), a = 1, p = l/2y/{3). 

In Figure 4 we plot the analytic solution (solid line) on a grid of 161 points. The dotted line is 
the numerical solution obtained with the Matlab routine ode45 with absolute and relative tolerance 
l.Oe — 4. The method produced this solution in 69 steps, and it was implemented with step-size 
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Figure 3: Error versus time in the SO (50) sparse case. 



control procedure. The dashed-dotted hne is the numerical solution obtained with the RK/MK 
method using SKC symmetric tecnhniques for the approximation of the exponential, with fixed 
step-size h. 

In Figure 5 we plot the error (along the y-axis) with respect to the numerical solution obtained 
with the Matlab routine ode45 to a logarithmic scale, versus the stepsize h = 1/2'° and k = 1, . . . ,5 
for the cases of the implementation of RK/MK with the expm function of Matlab (marked with +) 
and approximating exp to order four with a SKC technique (o). The line marked with * representes 
the error of the numerical solution given by the RK/MK method implemented with SKC technique 
for the approximation of the exponential, measured with respect to the numerical solution obtained 
by the same method with the use of the exact exponential (expm routine of Matlab). 

It is interesting to note in this case that substituting the exact exponential with suitable fourth- 
order approximant does not lead to a significant deterioration in the quality of the RK/MK method 
and the overall error does not change much. Note that in the present case the primary variable is 
a vector, rather than a matrix. In general, if the underlying ODE can be written in a vector form, 
i.e. as an action of a Lie group on R", we need to approximate ex-p{tB)v, where v G M", rather 
than the matrix exp(t_B). This leads to obvious savings in the SKC techniques, similarly, say, to the 
approach of rational functions. In particular, the cost of composing exponentials is 0(n^^, rather 
than C'(n'^), operations. 
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numerical solution and analytic soliUion 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 



Figure 4: The soliton originating in the KdV equation. 




Figure 5: RK/MK: global error at t = 5 with expm and SKC 
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A Appendix 

For completeness, we present a proof of Lemma 1. Note that a comprehensive treatment of this 
subject matter, inclusive of the non-symmetric case, has been presented in a different context by 
Chacon & Fomenko (1991). For our purposes, however, it is sufficient to derive the first term of the 
expansion. 

Lemma 1 The leading error term in the Strang splitting is 

s 

^ ^[Ci + • • • + G_i + iQ, [Ci + • • • + Ci-uQ]]. (4.6) 

1=2 

Proof Letting 

i^i(t) =e*^\ 

Fi{t) = c'^'^^Fi-ic'^'/^, l = 2,3,...,s, 
we can verify at once that Fg is precisely the Strang splitting. We assume that 

Flit) ^ cxp[t{Ci + --- + Ci) + ^Qit^ + 0{t^)], Z = 1, 2, . . . , s. 
We use the BCH formula: 

Flit) - exp(iiG) exp[t(C7i + • • • + G_i) + ^Qi-it^ + 0(<*)] exp(itCi) 

= exp{t(Ci + • • • + Q_i + iQ) + lt'-[Ci, Ci + • • • + 
+ - (Ci + • • • + [Q, Ci + • • • + 

+ ^t''Qi^i+Oit^)}eM¥Ci) 

= exp{f(Ci + ... + Ci) + lt'-[Ci,Ci + --- + 

+ - (Ci + • • • + C,_i), [G, Ci + • • • + C,_i]] 

+ + . . . + + iQ, G] + Ci + ■ • ■ + G] 

+ ^^^[^^ ^ . . . ^ [Ci + • • • + G„i + iQ, G]] + j^t^Q,-! + Oit^)}. 

However, 

[Cz, Ci + ■ • • + + [Ci + • • ■ + Q_i + iQ, Ci] = O, 
thereby annihilating the term, and 

^[iQ - (Ci + ■ ■ • + Q_i), [Cz,Ci + ■ ■ • + + j^[[Cz,Ci + • • • + Q_i],Ci] 
+ ^[Ci + • • • + [Ci + • • • + + ^Ci,Ci]] 
= + • • • + Q„i + iQ, [Ci + ■ ■ ■ + Q]]. 

Therefore 

Q; = [Ci + • • • + Q-i + iQ, [Ci + • • • + G]] + Q(_i. 
Since Qi = O, the expression (4.6) follows by summing the above formula for Z = 2, 3, . . . , s. □ 
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